pretemp    = pre ;
pre        = pre*36 ;
Postscore  = kappapostEB*36 ; 
DEPVAR_MeanSD = [mean(kappapost*36,'omitnan') std(kappapostEB*36,'omitnan')] %#ok
Deltakappa_MeanSD = [mean((kappapost-kappapre)*36,'omitnan') std((kappapostEB-kappapreEB)*36,'omitnan')] %#ok


[Ftemp,Xtemp]    = ecdf(AchievementData.Q(AchievementData.Q>=2)); Ftemp(2)=[]; Xtemp(2)=[] ;
tempbenchmark    = interp1(Ftemp,Xtemp,[(0.5-0.341); 0.5],'linear') ;
BenchmarkChangeQ = tempbenchmark(2)-tempbenchmark(1) %#ok
[Ftemp,Xtemp]    = ecdf(AchievementData.T(AchievementData.Q>=2)) ; Ftemp(2)=[]; Xtemp(2)=[] ; 
tempbenchmark    = interp1(Ftemp,Xtemp,[(0.5-0.341); 0.5],'linear') ;
BenchmarkChangeT = tempbenchmark(2)-tempbenchmark(1) %#ok


Tmean   = mean(T) ;
Tsig    = std(T) ;
T       = (T)/Tsig ; %%%%Here we normalize T for numerical stability
BenchmarkChangeTstd = (BenchmarkChangeT)/Tsig ;
Qmean   = mean(Q) ;
Qsig    = std(Q) ;
Q       = (Q)/Qsig ;
BenchmarkChangeQstd = (BenchmarkChangeQ)/Qsig ;
Premean = mean(pre) ;
Presig  = std(pre) ;
prestdz = (pre-Premean)/Presig ; %%%%We use standardized pre-test score in this section for numerical stability in Model 6
disp('******************************************************************************************************') 
disp('***********   HETEROSKEDASTICITY ROBUST ANALYSIS: SHORT-RUN HC PRODUCTION   **************************') 
disp('******************************************************************************************************') 
StDevDeltaScore  = std((kappapostEB-kappapreEB)*36,'omitnan') %#ok ; % std(Postscore,'omitnan') ; 
PostscoreMultiplier = StDevDeltaScore/std(Postscore,'omitnan') %#ok


ExpVars_SR5 = [ones(size(Postscore)), ... 1
               T, ... 2
               T.^2, ... 3
               Q, ... 4
               Q.^2, ...5
               T.*Q, ...6
               prestdz, ...7 %%%%We use standardized pre-test score for numerical stability in Model 6
               log(Y_EimputeShrunk), ...8
               log(Y_LimputeShrunk), ...9
               D2, ...10
               D3, ...11
               D2.*T, ...12
               D3.*T, ...13
               D2.*Q, ...14
               D3.*Q, ...15
               D2.*T.^2, ...16
               D3.*T.^2, ...17
               D2.*Q.^2, ...18
               D3.*Q.^2, ...19
               D2.*T.*Q, ...20
               D3.*T.*Q, ...21
               log(Y_EimputeShrunk).*T, ...22
               log(Y_LimputeShrunk).*T, ...23
               log(Y_EimputeShrunk).*Q, ...24
               log(Y_LimputeShrunk).*Q, ...25
               log(Y_EimputeShrunk).*T.^2, ...26
               log(Y_LimputeShrunk).*T.^2, ...27
               log(Y_EimputeShrunk).*Q.^2, ...28
               log(Y_LimputeShrunk).*Q.^2, ...29
               log(Y_EimputeShrunk).*T.*Q, ...30
               log(Y_LimputeShrunk).*T.*Q, ...31
               (prestdz).*T, ...32    %%%%We use standardized pre-test score for numerical stability in Model 6
               (prestdz).*T.^2, ...33 %%%%We use standardized pre-test score for numerical stability in Model 6
               (prestdz).*Q, ...34    %%%%We use standardized pre-test score for numerical stability in Model 6
               (prestdz).*Q.^2, ...35 %%%%We use standardized pre-test score for numerical stability in Model 6
               (prestdz).*T.*Q, ...36 %%%%We use standardized pre-test score for numerical stability in Model 6
               grade5, ...37
               grade5.*T, ...38
               grade5.*T.^2, ...39
               grade5.*Q, ...40
               grade5.*Q.^2, ...41
               grade5.*T.*Q, ... 42
               fem, ...43
               fem.*T, ...44
               fem.*T.^2, ...45
               fem.*Q, ...46
               fem.*Q.^2, ...47
               fem.*T.*Q, ...48
               blk, ...49
               blk.*T, ...50
               blk.*T.^2, ...51
               blk.*Q, ...52
               blk.*Q.^2, ...53
               blk.*T.*Q, ...54
               hsp, ...55
               hsp.*T, ...56
               hsp.*T.^2, ...57
               hsp.*Q, ...58
               hsp.*Q.^2, ...59
               hsp.*T.*Q, ...60
               SSEInc, ...61
               SSEInc.*T, ...62
               SSEInc.*T.^2, ...63
               SSEInc.*Q, ...64
               SSEInc.*Q.^2, ...65
               SSEInc.*T.*Q, ...66
               SSEIns, ...67
               SSEIns.*T, ...68
               SSEIns.*T.^2, ...69
               SSEIns.*Q, ...70
               SSEIns.*Q.^2, ...71
               SSEIns.*T.*Q, ...72
               NHelpers, ...73
               NHelpers.*T, ...74
               NHelpers.*T.^2, ...75
               NHelpers.*Q, ...76
               NHelpers.*Q.^2, ...77
               NHelpers.*T.*Q,... 78
               InvPar, ... 79
               InvPar.*T, ... 80
               InvPar.*T.^2, ... 81
               InvPar.*Q, ... 82
               InvPar.*Q.^2, ... 83
               InvPar.*T.*Q, ... 84
               BigFam, ... 85
               BigFam.*T, ... 86
               BigFam.*T.^2, ... 87
               BigFam.*Q, ... 88
               BigFam.*Q.^2, ... 89
               BigFam.*T.*Q, ... 90
               Middle, ... 91
               Middle.*T, ... 92
               Middle.*T.^2, ... 93
               Middle.*Q, ... 94
               Middle.*Q.^2, ... 95
               Middle.*T.*Q, ... 96
               Youngest,... 97
               Youngest.*T, ... 98
               Youngest.*T.^2, ... 99
               Youngest.*Q, ... 100
               Youngest.*Q.^2, ... 101
               Youngest.*T.*Q,... 102
               NoHmIntrnt, ... 103
               NoHmIntrnt.*T, ... 104
               NoHmIntrnt.*T.^2, ... 105
               NoHmIntrnt.*Q, ... 106
               NoHmIntrnt.*Q.^2, ... 107
               NoHmIntrnt.*T.*Q, ... 108
               MoblFrac, ... 109
               MoblFrac.*T, ... 110
               MoblFrac.*T.^2, ... 111
               MoblFrac.*Q, ... 112
               MoblFrac.*Q.^2, ... 113
               MoblFrac.*T.*Q, ... 114
               TabltFrac, ... 115
               TabltFrac.*T, ... 116
               TabltFrac.*T.^2, ... 117
               TabltFrac.*Q, ... 118
               TabltFrac.*Q.^2, ... 119
               TabltFrac.*T.*Q] ; % 120
%%
labels_SR5  = {'Const*';
               'T';'T^2'; ...
               'Q';'Q^2';'T*Q'; ...
               'PreScore*';'log(thetaE)';'log(thetaL)'; ...
               'D2';'D3'; 'D2*T'; 'D3*T'; 'D2*Q'; 'D3*Q'; 'D2*T^2'; 'D3*T^2'; 'D2*Q^2'; 'D3*Q^2'; 'D2*T*Q'; 'D3*T*Q'; ...
               'log(thetaE)*T'; 'log(thetaL)*T'; 'log(thetaE)*Q'; 'log(thetaL)*Q'; 'log(thetaE)*T^2'; 'log(thetaL)*T^2'; 'log(thetaE)*Q^2'; 'log(thetaL)*Q^2'; 'log(thetaE)*T*Q'; 'log(thetaL)*T*Q'; ...
               'pre*T';'pre*T^2';'pre*Q';'pre*Q^2';'pre*T*Q';'grade5';'grade5*T';'grade5*T^2';'grade5*Q';'grade5*Q^2';'grade5*T*Q';...
               'fem';'fem*T';'fem*T^2';'fem*Q';'fem*Q^2';'fem*T*Q';'blk';'blk*T';'blk*T^2';'blk*Q';'blk*Q^2';'blk*T*Q';'hsp';'hsp*T';'hsp*T^2';'hsp*Q';'hsp*Q^2';'hsp*T*Q';
               'SSEInc';'SSEInc*T';'SSEInc*T^2';'SSEInc*Q';'SSEInc*Q^2';'SSEInc*T*Q';'SSEIns';'SSEIns*T';'SSEIns*T^2';'SSEIns*Q';'SSEIns*Q^2';'SSEIns*T*Q';...
               '# HELPERS';'# HELPERS*T';'# HELPERS*T^2';'# HELPERS*Q';'# HELPERS*Q^2';'# HELPERS*T*Q';...
               'InvPar';'InvPar*T';'InvPar*T^2';'InvPar*Q';'InvPar*Q^2';'InvPar*T*Q';...
               'BigFam';'BigFam*T';'BigFam*T^2';'BigFam*Q';'BigFam*Q^2';'BigFam*T*Q';...
               'Middle';'Middle*T';'Middle*T^2';'Middle*Q';'Middle*Q^2';'Middle*T*Q';...
               'Youngest';'Youngest*T';'Youngest*T^2';'Youngest*Q';'Youngest*Q^2';'Youngest*T*Q';...
               'NoHmIntrnt';'NoHmIntrnt*T';'NoHmIntrnt*T^2';'NoHmIntrnt*Q';'NoHmIntrnt*Q^2';'NoHmIntrnt*T*Q';...
               'MoblFrac';'MoblFrac*T';'MoblFrac*T^2';'MoblFrac*Q';'MoblFrac*Q^2';'MoblFrac*T*Q';...
               'TabltFrac';'TabltFrac*T';'TabltFrac*T^2';'TabltFrac*Q';'TabltFrac*Q^2';'TabltFrac*T*Q'} ;
[alpha_SR5,se_SR5,tstat_SR5,pval_SR5,CI_SR5,out_SR5] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,[],[],0,[],[],[]); 
%%
StDevEffectATE    = nan(length(alpha_SR5)+2,1) ;
Delta0 = (alpha_SR5(1) + log(Y_EimputeShrunk)*alpha_SR5(8) + log(Y_LimputeShrunk)*alpha_SR5(9) ...
          + D2*alpha_SR5(10) + D3*alpha_SR5(11) + grade5*alpha_SR5(37) + fem*alpha_SR5(43) + blk*alpha_SR5(49) + hsp*alpha_SR5(55) ...
          + SSEInc*alpha_SR5(61) + SSEIns*alpha_SR5(67) + NHelpers*alpha_SR5(73) + InvPar*alpha_SR5(79) + BigFam*alpha_SR5(85) + Middle*alpha_SR5(91) + Youngest*alpha_SR5(97) ...
          + NoHmIntrnt*alpha_SR5(103) + MoblFrac*alpha_SR5(109) + TabltFrac*alpha_SR5(115) ) ;
StDevEffectATE(1) = std( Delta0 )/StDevDeltaScore ;  %%%%This is a mean SD effect of Delta0

dDeltadT = (alpha_SR5(2) + D2*alpha_SR5(12) + D3*alpha_SR5(13) + log(Y_EimputeShrunk)*alpha_SR5(22) + log(Y_LimputeShrunk)*alpha_SR5(23) ...
              +prestdz*alpha_SR5(32) + grade5*alpha_SR5(38) + fem*alpha_SR5(44) ...
              +blk*alpha_SR5(50) + hsp*alpha_SR5(56) + SSEInc*alpha_SR5(62) + SSEIns*alpha_SR5(68) + NHelpers*alpha_SR5(74) ...
              + InvPar*alpha_SR5(80) + BigFam*alpha_SR5(86) + Middle*alpha_SR5(92) + Youngest*alpha_SR5(98) ...
              + NoHmIntrnt*alpha_SR5(104) + MoblFrac*alpha_SR5(110) + TabltFrac*alpha_SR5(116) ) ...
            + 2*(alpha_SR5(3) + D2*alpha_SR5(16) + D3*alpha_SR5(17) + log(Y_EimputeShrunk)*alpha_SR5(26) + log(Y_LimputeShrunk)*alpha_SR5(27) ...
              + prestdz*alpha_SR5(33) + grade5*alpha_SR5(39) + fem*alpha_SR5(45) ...
              + blk*alpha_SR5(51) + hsp*alpha_SR5(57) + SSEInc*alpha_SR5(63) + SSEIns*alpha_SR5(69) + NHelpers*alpha_SR5(75) ...
              + InvPar*alpha_SR5(81) + BigFam*alpha_SR5(87) + Middle*alpha_SR5(93) + Youngest*alpha_SR5(99) ...
              + NoHmIntrnt*alpha_SR5(105) + MoblFrac*alpha_SR5(111) + TabltFrac*alpha_SR5(117) ).*T ...
            + Q.*(alpha_SR5(6) + D2*alpha_SR5(20) + D3*alpha_SR5(21) + log(Y_EimputeShrunk)*alpha_SR5(30) + log(Y_LimputeShrunk)*alpha_SR5(31) ...
              + prestdz*alpha_SR5(36) + grade5*alpha_SR5(42) + fem*alpha_SR5(48) ...
              + blk*alpha_SR5(54) + hsp*alpha_SR5(60) + SSEInc*alpha_SR5(66) + SSEIns*alpha_SR5(72) + NHelpers*alpha_SR5(78) ...
              + InvPar*alpha_SR5(84) + BigFam*alpha_SR5(90) + Middle*alpha_SR5(96) + Youngest*alpha_SR5(102) ...
              + NoHmIntrnt*alpha_SR5(108) + MoblFrac*alpha_SR5(114) + TabltFrac*alpha_SR5(120) ) ;
StDevEffectATE(2) = mean( BenchmarkChangeTstd*dDeltadT/StDevDeltaScore ) ;  %%%%This is a mean pseudo-SD effect of T

dDeltadQ = (alpha_SR5(4)+D2*alpha_SR5(14)+D3*alpha_SR5(15)+log(Y_EimputeShrunk)*alpha_SR5(24)+log(Y_LimputeShrunk)*alpha_SR5(25) ...
              +prestdz*alpha_SR5(34)+grade5*alpha_SR5(40)+fem*alpha_SR5(46) ...
              +blk*alpha_SR5(52)+hsp*alpha_SR5(58)+SSEInc*alpha_SR5(64)+SSEIns*alpha_SR5(70)+NHelpers*alpha_SR5(76) ...
              + InvPar*alpha_SR5(82) + BigFam*alpha_SR5(88) + Middle*alpha_SR5(94) + Youngest*alpha_SR5(100) ...
              + NoHmIntrnt*alpha_SR5(106) + MoblFrac*alpha_SR5(112) + TabltFrac*alpha_SR5(118) ) ...
            + 2*(alpha_SR5(5)+D2*alpha_SR5(18)+D3*alpha_SR5(19)+log(Y_EimputeShrunk)*alpha_SR5(28)+log(Y_LimputeShrunk)*alpha_SR5(29) ...
              +prestdz*alpha_SR5(35)+grade5*alpha_SR5(41)+fem*alpha_SR5(47) ...
              +blk*alpha_SR5(53)+hsp*alpha_SR5(59)+SSEInc*alpha_SR5(65)+SSEIns*alpha_SR5(71)+NHelpers*alpha_SR5(77) ...
              + InvPar*alpha_SR5(83) + BigFam*alpha_SR5(89) + Middle*alpha_SR5(95) + Youngest*alpha_SR5(101) ...
              + NoHmIntrnt*alpha_SR5(107) + MoblFrac*alpha_SR5(113) + TabltFrac*alpha_SR5(119) ).*Q ...
            + T.*(alpha_SR5(6)+D2*alpha_SR5(20)+D3*alpha_SR5(21)+log(Y_EimputeShrunk)*alpha_SR5(30)+log(Y_LimputeShrunk)*alpha_SR5(31) ...
              +prestdz*alpha_SR5(36)+grade5*alpha_SR5(42)+fem*alpha_SR5(48) ...
              +blk*alpha_SR5(54)+hsp*alpha_SR5(60)+SSEInc*alpha_SR5(66)+SSEIns*alpha_SR5(72)+NHelpers*alpha_SR5(78) ...
              + InvPar*alpha_SR5(84) + BigFam*alpha_SR5(90) + Middle*alpha_SR5(96) + Youngest*alpha_SR5(102) ...
              + NoHmIntrnt*alpha_SR5(108) + MoblFrac*alpha_SR5(114) + TabltFrac*alpha_SR5(120) ) ;
StDevEffectATE(4) = mean( BenchmarkChangeQstd*dDeltadQ/StDevDeltaScore ) ;  %%%%This is a standard deviation change in Q

StDevEffectATE(7) = mean( ( alpha_SR5(7) + T*alpha_SR5(32) + (T.^2)*alpha_SR5(33) + Q*alpha_SR5(34) + (Q.^2)*alpha_SR5(35) + T.*Q*alpha_SR5(36) - Presig )/StDevDeltaScore ) ;  %%%%This is the effect of a standard deviation change in Pre-Score (NOTE THAT PRE-SCORE IS STANDARDIZED, SO ITS STDEV IS 1).
StDevEffectATE(8) = mean( std(log(Y_EimputeShrunk))*( alpha_SR5(8) + T*alpha_SR5(22) + Q*alpha_SR5(24) + (T.^2)*alpha_SR5(26)+ (Q.^2)*alpha_SR5(28) + T.*Q*alpha_SR5(30) )/StDevDeltaScore ) ;  %%%%This is a standard deviation change in log(thetaE)
StDevEffectATE(9) = mean( std(log(Y_LimputeShrunk))*( alpha_SR5(9) + T*alpha_SR5(23) + Q*alpha_SR5(25) + (T.^2)*alpha_SR5(27)+ (Q.^2)*alpha_SR5(29) + T.*Q*alpha_SR5(31) )/StDevDeltaScore ) ;  %%%%This is a standard deviation change in log(thetaL)
StDevEffectATE(10) = mean( ( alpha_SR5(10) + T*alpha_SR5(12) + Q*alpha_SR5(14) + (T.^2)*alpha_SR5(16) + (Q.^2)*alpha_SR5(18) + T.*Q*alpha_SR5(20) )/StDevDeltaScore ) ;  %%%%This is a change from District 1 to District 2
StDevEffectATE(11) = mean( ( alpha_SR5(11) + T*alpha_SR5(13) + Q*alpha_SR5(15) + (T.^2)*alpha_SR5(17) + (Q.^2)*alpha_SR5(19) + T.*Q*alpha_SR5(21) )/StDevDeltaScore ) ;  %%%%This is a change from District 1 to District 3
StDevEffectATE(37) = mean( ( alpha_SR5(37) + T*alpha_SR5(38) + (T.^2)*alpha_SR5(39) + Q*alpha_SR5(40) + (Q.^2)*alpha_SR5(41) + T.*Q*alpha_SR5(42) )/StDevDeltaScore ) ;  %%%%This is a change from grade 6 to grade 5
StDevEffectATE(43) = mean( ( alpha_SR5(43) + T*alpha_SR5(44) + (T.^2)*alpha_SR5(45) + Q*alpha_SR5(46) + (Q.^2)*alpha_SR5(47) + T.*Q*alpha_SR5(48) )/StDevDeltaScore ) ;  %%%%This is a change from male to female
StDevEffectATE(49) = mean( ( alpha_SR5(49) + T*alpha_SR5(50) + (T.^2)*alpha_SR5(51) + Q*alpha_SR5(52) + (Q.^2)*alpha_SR5(53) + T.*Q*alpha_SR5(54) )/StDevDeltaScore ) ;  %%%%This is a change from wao to black
StDevEffectATE(55) = mean( ( alpha_SR5(55) + T*alpha_SR5(56) + (T.^2)*alpha_SR5(57) + Q*alpha_SR5(58) + (Q.^2)*alpha_SR5(59) + T.*Q*alpha_SR5(60) )/StDevDeltaScore ) ;  %%%%This is a change from wao to hispanic
StDevEffectATE(61) = mean( std(SSEInc)*( alpha_SR5(61) + T*alpha_SR5(62) + (T.^2)*alpha_SR5(63) + Q*alpha_SR5(64) + (Q.^2)*alpha_SR5(65) + T.*Q*alpha_SR5(66) )/StDevDeltaScore ) ;  %%%%This is a StDev change in SSEInc
StDevEffectATE(67) = mean( std(SSEIns)*( alpha_SR5(67) + T*alpha_SR5(68) + (T.^2)*alpha_SR5(69) + Q*alpha_SR5(70) + (Q.^2)*alpha_SR5(71) + T.*Q*alpha_SR5(72) )/StDevDeltaScore ) ;  %%%%This is a StDev change in SSEIns
StDevEffectATE(73) = mean( std(NHelpers)*( alpha_SR5(73) + T*alpha_SR5(74) + (T.^2)*alpha_SR5(75) + Q*alpha_SR5(76) + (Q.^2)*alpha_SR5(77) + T.*Q*alpha_SR5(78) )/StDevDeltaScore ) ;  %%%%This is a StDev change in NHelpers

OtherEffect    = nan(length(alpha_SR5)+2,1) ;
%%%%This is a mean total derivative of kappa with respect to Q, for a
%%%%change of BenchmarkChangeQstd, as a fraction of StDevDeltaScore
OtherEffect(4) = mean( ( dDeltadQ + dDeltadT.*( (EPointEst.tau0*Y_Eimpute.*((Q*Qsig+1).^(-EPointEst.phi)))/Tsig ) )*BenchmarkChangeQstd/StDevDeltaScore ) ;  %%%%This is a mean total derivative effect of Q (in original DeltaS units)
OtherEffect(7) = mean( ( T*alpha_SR5(32) + (T.^2)*alpha_SR5(33) + Q*alpha_SR5(34) + (Q.^2)*alpha_SR5(35) + T.*Q*alpha_SR5(36) )/StDevDeltaScore ) ;  %%%%This is the effect of a standard deviation change in Pre-Score SLOPE EFFECTS ONLY (NOTE THAT PRE-SCORE IS STANDARDIZED, SO ITS STDEV IS 1).
OtherEffect(8) = mean( std(log(Y_EimputeShrunk))*( T*alpha_SR5(22) + Q*alpha_SR5(24) + (T.^2)*alpha_SR5(26)+ (Q.^2)*alpha_SR5(28) + T.*Q*alpha_SR5(30) )/StDevDeltaScore ) ;  %%%%This is a standard deviation change in log(thetaE) SLOPE EFFECTS ONLY
OtherEffect(9) = mean( std(log(Y_LimputeShrunk))*( T*alpha_SR5(23) + Q*alpha_SR5(25) + (T.^2)*alpha_SR5(27)+ (Q.^2)*alpha_SR5(29) + T.*Q*alpha_SR5(31) )/StDevDeltaScore ) ;  %%%%This is a standard deviation change in log(thetaL) SLOPE EFFECTS ONLY
OtherEffect(10) = mean( ( T*alpha_SR5(12) + Q*alpha_SR5(14) + (T.^2)*alpha_SR5(16) + (Q.^2)*alpha_SR5(18) + T.*Q*alpha_SR5(20) )/StDevDeltaScore ) ;  %%%%This is a change from District 1 to District 2 SLOPE EFFECTS ONLY
OtherEffect(11) = mean( ( T*alpha_SR5(13) + Q*alpha_SR5(15) + (T.^2)*alpha_SR5(17) + (Q.^2)*alpha_SR5(19) + T.*Q*alpha_SR5(21) )/StDevDeltaScore ) ;  %%%%This is a change from District 1 to District 3 SLOPE EFFECTS ONLY
OtherEffect(12) = std( ( T*alpha_SR5(12) + Q*alpha_SR5(14) + (T.^2)*alpha_SR5(16) + (Q.^2)*alpha_SR5(18) + T.*Q*alpha_SR5(20) ) ) ;  %%%%This is a change from District 1 to District 2 SLOPE EFFECTS ONLY
OtherEffect(13) = std( ( T*alpha_SR5(13) + Q*alpha_SR5(15) + (T.^2)*alpha_SR5(17) + (Q.^2)*alpha_SR5(19) + T.*Q*alpha_SR5(21) ) ) ;  %%%%This is a change from District 1 to District 3 SLOPE EFFECTS ONLY




disp('-----------------------------------------------------------------------------------')
disp('-------------------------  SHORT-RUN PROGRESS: MODEL 5  ---------------------------')
disp('-----------------------------------------------------------------------------------')
baselinetemp = sort( alpha_SR5(1)+prestdz*alpha_SR5(7)+log(Y_EimputeShrunk)*alpha_SR5(8)+log(Y_LimputeShrunk)*alpha_SR5(9) ...
                            +D2*alpha_SR5(10)+D3*alpha_SR5(11)+grade5*alpha_SR5(37)+fem*alpha_SR5(43)+blk*alpha_SR5(49)+hsp*alpha_SR5(55) ...
                            +SSEInc*alpha_SR5(61)+SSEIns*alpha_SR5(67)+NHelpers*alpha_SR5(73) + InvPar*alpha_SR5(79) + BigFam*alpha_SR5(85) + Middle*alpha_SR5(91) + Youngest*alpha_SR5(97) ...
                            + NoHmIntrnt*alpha_SR5(103) + MoblFrac*alpha_SR5(109) + TabltFrac*alpha_SR5(115) - pre ) ;
BASELINE_SR5 = [mean(baselinetemp) median(baselinetemp) std(baselinetemp) baselinetemp(round(length(baselinetemp)*0.1)) baselinetemp(round(length(baselinetemp)*0.9))] %#ok
out_SR5.results.StDevEffectATE = StDevEffectATE;
out_SR5.results.OtherEffect = OtherEffect;
WhiteTest  = out_SR5.White %#ok
Method     = out_SR5.method %#ok
disp(out_SR5.results(:,:)) 
%%
%%%%Test for joint significance of all intercept terms:           
R = zeros(20,length(ExpVars_SR5(1,:))) ; R(1,1)=1; R(2,7)=1; R(3,8)=1; R(4,9)=1; R(5,10)=1; R(6,11)=1;... 
                                         R(7,37)=1; R(8,43)=1; R(9,49)=1; R(10,55)=1; R(11,61)=1; R(12,67)=1;... 
                                         R(13,73)=1; R(14,79)=1; R(15,85)=1; R(16,91)=1; R(17,97)=1; R(18,103)=1; ...
                                         R(19,109)=1; R(20,115)=1; r = zeros(20,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
InterceptTest = temp.Wald %#ok
%%%%Test for joint significance of all NON-S^{pre} intercept terms:           
R = zeros(19,length(ExpVars_SR5(1,:))) ; R(1,1)=1; R(2,8)=1; R(3,9)=1; R(4,10)=1; R(5,11)=1; R(6,37)=1; ...
                                         R(7,43)=1; R(8,49)=1; R(9,55)=1; R(10,61)=1; R(11,67)=1; R(12,73)=1; ...
                                         R(13,79)=1; R(14,85)=1; R(15,91)=1; R(16,97)=1; R(17,103)=1; R(18,109)=1; ...
                                         R(19,115)=1; r = zeros(19,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
InterceptTest_NonSpre = temp.Wald %#ok
%%%%Test for joint significance of all T terms:           
R = zeros(60,length(ExpVars_SR5(1,:))) ; R(1,2)=1; R(2,3)=1; R(3,6)=1; R(4,12)=1; R(5,13)=1; R(6,16)=1; ...
                                         R(7,17)=1; R(8,20)=1; R(9,21)=1; R(10,22)=1; R(11,23)=1; R(12,26)=1; 
                                         R(13,27)=1; R(14,30)=1; R(15,31)=1; R(16,32)=1; R(17,33)=1; R(18,36)=1;
                                         R(19,38)=1; R(20,39)=1; R(21,42)=1; R(22,44)=1; R(23,45)=1; R(24,48)=1;
                                         R(25,50)=1; R(26,51)=1; R(27,54)=1; R(28,56)=1; R(29,57)=1; R(30,60)=1;
                                         R(31,62)=1; R(32,63)=1; R(33,66)=1; R(34,68)=1; R(35,69)=1; R(36,72)=1;
                                         R(37,74)=1; R(38,75)=1; R(39,78)=1; R(40,80)=1; R(41,81)=1; R(42,84)=1;
                                         R(43,86)=1; R(43,87)=1; R(44,90)=1; R(45,92)=1; R(46,93)=1; R(47,96)=1; 
                                         R(43,98)=1; R(43,99)=1; R(48,102)=1; R(49,104)=1; R(50,105)=1; R(51,108)=1; 
                                         R(52,110)=1; R(53,111)=1; R(54,114)=1; R(55,116)=1; R(56,117)=1; R(43,120)=1; 
                                         r = zeros(60,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
TTest = temp.Wald %#ok
%%%%Test for joint significance of all Q terms:           
R = zeros(60,length(ExpVars_SR5(1,:))) ; R(1,4)=1; R(2,5)=1; R(3,6)=1; R(4,14)=1; R(5,15)=1; R(6,18)=1; 
                                         R(7,19)=1; R(8,20)=1; R(9,21)=1; R(10,24)=1; R(11,25)=1; R(12,28)=1; 
                                         R(13,29)=1; R(14,30)=1; R(15,31)=1; R(16,34)=1; R(17,35)=1; R(18,36)=1;
                                         R(19,40)=1; R(20,41)=1; R(21,42)=1; R(22,46)=1; R(23,47)=1; R(24,48)=1;
                                         R(25,52)=1; R(26,53)=1; R(27,54)=1; R(28,58)=1; R(29,59)=1; R(30,60)=1; 
                                         R(31,64)=1; R(32,65)=1; R(33,66)=1; R(34,70)=1; R(35,71)=1; R(36,72)=1; 
                                         R(37,76)=1; R(38,77)=1; R(39,78)=1; R(40,82)=1; R(41,83)=1; R(42,84)=1;
                                         R(43,88)=1; R(44,89)=1; R(45,90)=1; R(46,94)=1; R(47,95)=1; R(48,96)=1;
                                         R(49,100)=1; R(50,101)=1; R(51,102)=1; R(52,106)=1; R(53,107)=1; R(54,108)=1; 
                                         R(55,112)=1; R(56,113)=1; R(55,114)=1; R(56,118)=1; R(55,119)=1; R(56,120)=1; 
                                         r = zeros(60,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
QTest = temp.Wald %#ok   
%%%%Test for joint significance of all District 2 terms:           
R = zeros(6,length(ExpVars_SR5(1,:))) ; R(1,10)=1; R(2,12)=1; R(3,14)=1; R(4,16)=1; R(5,18)=1; R(6,20)=1; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
Dist2Test = temp.Wald %#ok
%%%%Test for joint significance of all District 2 slope interactions:           
R = zeros(5,length(ExpVars_SR5(1,:))) ; R(1,12)=1; R(2,14)=1; R(3,16)=1; R(4,18)=1; R(5,20)=1; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
Dist2InteractionsTest = temp.Wald %#ok 
%%%%Test for joint significance of all District 3 terms:           
R = zeros(6,length(ExpVars_SR5(1,:))) ; R(1,11)=1; R(2,13)=1; R(3,15)=1; R(4,17)=1; R(5,19)=1; R(6,21)=1; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
Dist3Test = temp.Wald %#ok
%%%%Test for joint significance of all District 3 slope interactions:           
R = zeros(5,length(ExpVars_SR5(1,:))) ; R(1,13)=1; R(2,15)=1; R(3,17)=1; R(4,19)=1; R(5,21)=1; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
Dist3InteractionsTest = temp.Wald %#ok 
%%%%Test for joint significance of all School District terms:           
R = zeros(12,length(ExpVars_SR5(1,:))) ; for ii=1:12; R(ii,9+ii) = 1 ; end; r = zeros(12,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
DistTest = temp.Wald %#ok
%%%%Test for joint significance of all School District slope interactions:           
R = zeros(10,length(ExpVars_SR5(1,:))) ; for ii=1:10; R(ii,11+ii) = 1 ; end; r = zeros(10,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
DistInteractionsTest = temp.Wald %#ok 
%%%%Test for joint significance of all School District intercept terms:           
R = zeros(2,length(ExpVars_SR5(1,:))) ; R(1,10) = 1 ; R(2,11) = 1 ;  r = zeros(2,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
DistInterceptTest = temp.Wald %#ok       
%%%%Test for joint significance of all ThetaE terms: 
R = zeros(6,length(ExpVars_SR5(1,:))) ; R(1,8)=1 ; R(2,22)=1 ; R(3,24)=1 ; R(4,26)=1 ; R(5,28)=1 ; R(6,30)=1 ; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
ThetaETest = temp.Wald %#ok           
%%%%Test for joint significance of all ThetaE interaction terms: 
R = zeros(5,length(ExpVars_SR5(1,:))) ; R(1,22)=1 ; R(2,24)=1 ; R(3,26)=1 ; R(4,28)=1 ; R(5,30)=1 ; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
ThetaEInteractionsTest = temp.Wald %#ok 
%%%%Test for joint significance of all ThetaL terms:          
R = zeros(6,length(ExpVars_SR5(1,:))) ; R(1,9)=1 ; R(2,23)=1 ; R(3,25)=1 ; R(4,27)=1 ; R(5,29)=1 ; R(6,31)=1 ; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
ThetaLTest = temp.Wald %#ok 
%%%%Test for joint significance of all ThetaL interaction terms:          
R = zeros(5,length(ExpVars_SR5(1,:))) ; R(1,23)=1 ; R(2,25)=1 ; R(3,27)=1 ; R(4,29)=1 ; R(5,31)=1 ; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
ThetaLInteractionsTest = temp.Wald %#ok 
%%%%Test for joint significance of all Student Type slope interactions:          
R = zeros(10,length(ExpVars_SR5(1,:))) ; for ii=1:10; R(ii,21+ii) = 1 ; end; r = zeros(10,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
TypeInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all pre-test terms:           
R = zeros(6,length(ExpVars_SR5(1,:))) ; R(1,7)=1; for ii=1:5; R(ii+1,31+ii)=1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
PreScoreTest = temp.Wald %#ok
%%%%Test for joint significance of all pre-test slope interactions:           
R = zeros(5,length(ExpVars_SR5(1,:))) ; for ii=1:5; R(ii,31+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
PreScoreInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all age-cohort terms:           
R = zeros(6,length(ExpVars_SR5(1,:))) ; for ii=1:6; R(ii,36+ii) = 1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
Grade5Test = temp.Wald %#ok
%%%%Test for joint significance of all age-cohort slope interactions:           
R = zeros(5,length(ExpVars_SR5(1,:))) ; for ii=1:5; R(ii,37+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
Grade5InteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all fem terms:           
R = zeros(6,length(ExpVars_SR5(1,:))) ; for ii=1:6; R(ii,42+ii) = 1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
FemTest = temp.Wald %#ok
%%%%Test for joint significance of all fem slope interactions:           
R = zeros(5,length(ExpVars_SR5(1,:))) ; for ii=1:5; R(ii,43+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
FemInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all blk terms:           
R = zeros(6,length(ExpVars_SR5(1,:))) ; for ii=1:6; R(ii,48+ii) = 1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
BlkTest = temp.Wald %#ok
%%%%Test for joint significance of all blk slope interactions:           
R = zeros(5,length(ExpVars_SR5(1,:))) ; for ii=1:5; R(ii,49+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
BlkInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all hsp terms:           
R = zeros(6,length(ExpVars_SR5(1,:))) ; for ii=1:6; R(ii,54+ii) = 1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
HspTest = temp.Wald %#ok
%%%%Test for joint significance of all hsp slope interactions:           
R = zeros(5,length(ExpVars_SR5(1,:))) ; for ii=1:5; R(ii,55+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
HspInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all SSEInc terms:           
R = zeros(6,length(ExpVars_SR5(1,:))) ; for ii=1:6; R(ii,60+ii) = 1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
SSEIncTest = temp.Wald %#ok
%%%%Test for joint significance of all SSEInc slope interactions:           
R = zeros(5,length(ExpVars_SR5(1,:))) ; for ii=1:5; R(ii,61+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
SSEIncInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all SSEIns terms:           
R = zeros(6,length(ExpVars_SR5(1,:))) ; for ii=1:6; R(ii,66+ii) = 1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
SSEInsTest = temp.Wald %#ok
%%%%Test for joint significance of all SSEIns slope interactions:           
R = zeros(5,length(ExpVars_SR5(1,:))) ; for ii=1:5; R(ii,67+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
SSEInsInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all SSE terms:           
R = zeros(12,length(ExpVars_SR5(1,:))) ; for ii=1:12; R(ii,60+ii) = 1 ; end; r = zeros(12,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
SSETest = temp.Wald %#ok
%%%%Test for joint significance of all NHelper terms:           
R = zeros(6,length(ExpVars_SR5(1,:))) ; for ii=1:6; R(ii,72+ii) = 1 ; end; r = zeros(6,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
NHelperTest = temp.Wald %#ok
%%%%Test for joint significance of all NHelper slope interactions:           
R = zeros(5,length(ExpVars_SR5(1,:))) ; for ii=1:5; R(ii,73+ii) = 1 ; end; r = zeros(5,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
NHelperInteractionsTest = temp.Wald %#ok
%%%%Test for joint significance of all Parent Survey Terms:           
R = zeros(24,length(ExpVars_SR5(1,:))) ; for ii=1:24; R(ii,78+ii) = 1 ; end; r = zeros(24,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
ParSurveyTest = temp.Wald %#ok
%%%%Test for joint significance of all Home Connectivity Terms:           
R = zeros(18,length(ExpVars_SR5(1,:))) ; for ii=1:18; R(ii,102+ii) = 1 ; end; r = zeros(18,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
HomeConnectivityTest = temp.Wald %#ok
%%%%Test for joint significance of various T*Q Terms:           
R = zeros(11,length(ExpVars_SR5(1,:))) ; R(5,66) = 1 ; R(6,72) = 1 ; R(7,78) = 1 ; R(8,104) = 1 ; ...
                                         R(5,84) = 1 ; R(6,90) = 1 ; R(7,96) = 1 ; R(8,102) = 1 ; ...
                                         R(9,108) = 1 ; R(10,114) = 1 ; R(10,120) = 1 ; 
                                         r = zeros(11,1) ;
[~,~,~,~,~,temp] = regressHet(Postscore,ExpVars_SR5,'FGLSrobust',labels_SR5,R,r,0,[],[],[]); 
TtimesQTest = temp.Wald %#ok

                                








%%
Q = AchievementData.Q ; Q = Q(useindx) ;  %%%%THIS WILL BE USED AS A REGRESSOR IN THE HC PRODUCTION TECHNOLOGY EQUATIONS
T = AchievementData.T ; T = T(useindx) ;  %%%%THIS WILL BE USED AS A REGRESSOR IN THE HC PRODUCTION TECHNOLOGY EQUATIONS


pre = pretemp ;
diary off